Introduction

Thomas Wiecki created a very interesting blog post on the Markov chain Monte Carlo (MCMC) Metropolis sampling algorithm . Here I attempt to work through it and expand on it where I feel like.

The issue being addressed relates to the Bayes formula: $P(\theta|x) = \frac{P(x|\theta) P(\theta)}{P(x)}$

Here we are interested in the probability of our model parameters $\theta$ given the data $x$. Apparently calculating $P(x)$ (aka. the evidence (the evidence that $x$ was generated by this model)) is the real ballache requiring integration: $P(x) = \int_\Theta P(x, \theta) \, \mathrm{d}\theta$ So, it is very hard to solve for, and even harded to sample from. But using MCMC it is possible to construct a Markov chain to match the posterior distribution.

Goal

The gola for the explorations below is to estimate/infer the posterior of the mean $\mu$ assuming known standard deviation $(\sigma) = 1$.

Load libraries and set up environment



In [18]:

    
%matplotlib inline

import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from scipy.stats import norm

Generate data for experimentation

Data set

100 points from a normal centered around zero



In [19]:

    
np.random.seed(3141)
data = np.random.randn(20)



In [20]:

    
ax = plt.subplot()
sns.distplot(data, kde=False, ax=ax)
_ = ax.set(title='Histogram of observed data', xlabel='x', ylabel='# observations');

Model definition

This problem's model is:

$\mu \sim \text{Normal}(0, 1)\\ x|\mu \sim \text{Normal}(x; \mu, 1)$

Caclulate the the parameters of the posterior



In [21]:

    
def calc_posterior_analytical(data, x, mu_0, sigma_0):
    sigma = 1.
    n = len(data)
    mu_post = (mu_0 / sigma_0**2 + data.sum() / sigma**2) / (1. / sigma_0**2 + n / sigma**2)
    sigma_post = (1. / sigma_0**2 + n / sigma**2)**-1
    return norm(mu_post, np.sqrt(sigma_post)).pdf(x)

Show the probability of $\mu$'s values after having seen the data, taking prior information into account



In [22]:

    
ax = plt.subplot()
x = np.linspace(-1, 1, 80)
posterior_analytical = calc_posterior_analytical(data, x, 0., 1.)
ax.plot(x, posterior_analytical)
ax.set(xlabel='mu', ylabel='belief', title='Analytical posterior');
sns.despine()



In [ ]: